function grad = gradLL2_2(alpha,choiceFull,dPds,dPdv,dvdtheta,dzetadalpha,...
    dzetads,dzetadv,H,K,LL,phi,ppFull,tau,tcidx,P,pi,zeta)
% Calculate gradient of log likelihood (LL)

% Initialize parameters and data matrices----------------------------------
obs = size(choiceFull(:,:,1),1); % Number of observations
M = H - 1 + tau - 1 + 1; % Number of variables in theta = [mu eta chi]
dLdtheta = zeros(obs,M);
dLdpi = zeros(obs,tau-1);
dLds = zeros(obs,1);
dLdalpha = zeros(obs,1);


for idx1 = 1:obs
    
    % Collect choice and travel cost data for each observation
    tcCombo = tcidx(idx1);
    prefPts08 = ppFull{1}(idx1) + 1;
    prefPts09 = min(ppFull{2}(idx1) + 1,K);
    hunt08 = choiceFull(idx1,1) + 1;
    hunt09 = choiceFull(idx1,2) + 1;
    phi08 = phi(hunt08,prefPts08);
    if prefPts09 == 1 && hunt08 < 24
        G = phi08;
    elseif hunt08 == 24
        G = 1;
    elseif prefPts09 > 1 && hunt08 < 24
        G = (1 - phi08);
    end
          
    
    dLdalpha_num = zeros(1,tau);
    dLds_num = zeros(1,tau);
    dLdpi_num = zeros(1,tau-1);
    dLdtheta_num = zeros(tau,M);
    for idx2 = 1:tau
        
        % Initialize variables
        P09 = P(hunt09,prefPts09,tcCombo,idx2);
        P08 = P(hunt08,prefPts08,tcCombo,idx2);
        P09_3 = P(hunt09,prefPts09,tcCombo,3);
        P08_3 = P(hunt08,prefPts08,tcCombo,3);
        zeta08 = zeta(prefPts08,1,tcCombo,idx2);
        zeta08_3 = zeta(prefPts08,1,tcCombo,3);
        dzetadalpha_ = dzetadalpha(prefPts08,1,tcCombo,idx2);
        dPds09 = dPds(hunt09,prefPts09,tcCombo,idx2);
        dPds08 = dPds(hunt08,prefPts08,tcCombo,idx2);
        dzetads_ = dzetads(prefPts08,1,tcCombo,idx2);
        dPdv09 = dPdv(hunt09,:,prefPts09,tcCombo,idx2);
        dPdv08 = dPdv(hunt08,:,prefPts08,tcCombo,idx2);
        dzetadv_ = dzetadv(prefPts08,:,tcCombo,idx2);
        dvdtheta_ = dvdtheta(:,:,tcCombo,idx2);
        
        
        % Calculate dLL/dalpha
        dP1dalpha = P09*G*P08*(dzetadalpha_*(1 - alpha)/(1 + alpha)...
            - zeta08*2/(1 + alpha)^2);
        dP2dalpha = P08*(dzetadalpha_*alpha/(1 + alpha) + zeta08/(1 + alpha)^2);
        dP3dalpha = P09/(1 + alpha)^2;
        
        dLdalpha_num(1,idx2) = (dP1dalpha + (hunt09==24)*dP2dalpha...
            + (hunt08==24&prefPts09==1)*dP3dalpha)*pi(idx2);
        
        
        % Calculate dLL/ds
        dP1ds = G*(1 - alpha)*pi(idx2)/(1 + alpha)...
            *(dPds09*P08*zeta08 + P09*dPds08*zeta08 + P09*P08*dzetads_);
        dP2ds = alpha*pi(idx2)/(1 + alpha)*(dPds08*zeta08 + P08*dzetads_);
        dP3ds = dPds09*alpha/(1+alpha)*pi(idx2);
    
        dLds_num(1,idx2) = dP1ds + (hunt09==24)*dP2ds ...
            + (hunt08==24&prefPts09==1)*dP3ds;
        
        
        % Calculate dLL/dpi

        if idx2 < 3 

            dP1dpi = G*(1 - alpha)/(1 + alpha)...
                *(P09*P08*zeta08 - P09_3*P08_3*zeta08_3);
            dP2dpi = alpha/(1 + alpha)*(P08*zeta08 - P08_3*zeta08_3);
            dP3dpi = alpha/(1 + alpha)*(P09 - P09_3);

            dLdpi_num(1,idx2) = dP1dpi + (hunt09==24)*dP2dpi ...
                + (hunt08==24&prefPts09==1)*dP3dpi;
            
        end
        
        
        % Calculate dLL/dtheta, theta = [mu eta chi]
        dP1dtheta = G*(1 - alpha)*pi(idx2)/(1 + alpha)...
            *(dPdv09*dvdtheta_*P08*zeta08 + P09*dPdv08*dvdtheta_*zeta08...
            + P09*P08*dzetadv_*dvdtheta_);
        dP2dtheta = alpha*pi(idx2)/(1 + alpha)*(dPdv08*dvdtheta_*zeta08...
            + P08*dzetadv_*dvdtheta_);
        dP3dtheta = pi(idx2)*alpha/(1 + alpha)*dPdv09*dvdtheta_;
        
        dLdtheta_num(idx2,:) = dP1dtheta + (hunt09==24)*dP2dtheta ...
            + (hunt08==24&prefPts09==1)*dP3dtheta;
        
    end
    
    dLdalpha(idx1) = sum(dLdalpha_num)/exp(LL(idx1));
    dLds(idx1) = sum(dLds_num)/exp(LL(idx1));    
    dLdpi(idx1,:) = dLdpi_num./exp(LL(idx1));
    dLdtheta(idx1,:) = sum(dLdtheta_num)./exp(LL(idx1));
    
end

grad = -sum([dLdtheta(:,1:3) dLdpi dLdtheta(:,4:end) dLdalpha dLds]);












